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The quasi-equilibrium approximation is employed as a systematic tool for solving the problem of deriving 
constitutive equations from kinetic models of liquid-crystalline polymers. It is demonstrated how kinetic models 
of liquid-crystalline polymers can be approximated in a systematic way, how canonical distribution functions 
can be derived from the maximum entropy principle and how constitutive equations are derived therefrom. 
The numerical implementation of the constitutive equations based on the intrinsic dual structure of the quasi- 
equilibrium manifold thus derived is developed and illustrated for particular examples. Finally, a measure of 
the accuracy of the quasi-equilibrium approximation is proposed that can be implemented into the numerical 
integration of the constitutive equation. 
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to ' I- INTRODUCTION 

+-* . 

This paper continues a systematic approach to the derivation and numerical implementation of constitutive equations for 
complex fluids, initiated in Ref. In this approach, the systematic application of the quasi-equilibrium approximation to 

' kinetic models is proposed, including the derivation of canonical distribution functions (CDF) and their use to obtain closed 
form constitutive equations. Also the natural numerical implementation of the resulting constitutive equations based on the dual 
structure of the quasi-equilibria is discussed together with a measure of the accuracy of the quasi-equilibrium approximation. 
While in Ref. [Jj]] dilute solutions of flexible polymers were considered, the present paper extends this study to kinetic models of 
liquid-crystalline polymers. Peculiarities due to the mean-field nature of models for liquid crystalline polymers are discussed. 

| A more detailed presentation of the subject can be found in [EJ. 

. Liquid-crystalline polymers have attracted considerable attention due to their capability of forming a highly oriented phase 
Ml. In polymer processing, for example, it is very important to understand the interplay between the the tendency towards 
orientational ordering and the orientational effects due to the flow field [Q, |5j| . Kinetic models of liquid-crystalline (rigid rodlike) 
polymers were introduced by Hess and Doi |^, ^[ ||] and are now used by many authors. Kinetic models of the dynamics of liquid- 
crystalline polymers are able to reproduce qualitatively most of the experimental results for the steady state in homogeneous flows 
Predictions of the kinetic models for transient viscoelastic phenomena are in general less reliable jlO|]. The main limitation 
of the kinetic model responsible for this failure is probably the neglect of defects of the orientational ordering. Recently, some 
work has been done in order to implement spatial variations of the orientational order into the kinetic model flu] , |l2| ] . For clarity, 
we restrict ourselves here to spatially homogeneous systems, and the corresponding kinetic models. The present approach can 
£h , be easily extended to include also spatially inhomogeneous systems. 

Numerical implementation of kinetic models of liquid-crystalline polymers in direct numerical flow calculations is in general 
very computationally intensive. However, kinetic models of polymer dynamics may serve as a starting point for the derivation 
of constitutive equations. In general, this derivation is not straightforward but requires approximations to the underlying kinetic 
model. The need for so-called closure approximations occurs in many branches of statistical physics and several suggestions 
for such approximations have been proposed in the literature (see e.g. [ pj| ] and references therein). A particularly successful 
approach for the derivation of constitutive equation is the so-called method of invariant manifold Jl4h. This method identifies 
relevant manifolds for the reduced description from which the constitutive equations can be derived. For example, it was found 
in Ref. [|l5[ with the help of this method that the universal limit in the dynamics of dilute polymer solutions is given by the 
(revised) Oldroyd 8-constant model. 

In Ref. and in the present work, the quasi-equilibrium approximation is employed in order to derive a manifold formed by 
a set of canonical distribution functions (CDF). The quasi-equilibrium approximation is used successfullyalso in other branches 
of statistical physics, like for example chemical reaction kinetics, plasma and ferrofluids (see also Refs. jl^ , [l6|]). This approxi- 
mation to the dynamics shows several desirable features like conservation of the dissipative nature of the dynamics, conservation 
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of positive-definiteness of distribution functions and a thermodynamic structure including dual variables. Improvements on this 
approximation can be found by either including more variables into the QEA, by constructing improved manifolds according to 
the method of invariant manifold proposed in Ref. [ |l4| ] or by using a combined micro-macro simulation technique proposed in 
Ref. [[171], Here, we show how the quasi-equilibrium approximation can be used in order to obtain closed constitutive equations 
for liquid-crystalline polymers. We also present an algorithm for the numerical implementation of the constitutive equations. In 
addition, a measure for the accuracy of the approximation is suggested that does not require solutions to the kinetic model and 
can therefore be used while solving the constitutive equation. 

Recently, the idea of using a set of canonical distribution functions (CDF), which are postulated and not based on the quasi- 
equilibrium approximation, for obtaining closures to kinetic models of dilute solutions of flexible polymers has also been used 
in Refs. |J[ 

This paper is organized as follows. Kinetic models of the dynamics of liquid-crystalline polymers (LCPs) are reviewed in 
Sec. II. Canonical distribution functions and their use in the description of the dynamics of LCPs is described in Sec. |l| Closed 
form constitutive equations are derived using CDF. Special emphasize is paid to so-called alignment tensor models. In Sec. [v|, a 
measure for the accuracy of the approximate description of polymer dynamics with CDF is p roposed. In Sec. [VlJ some numerical 
results in steady shear flow are presented. Finally, some conclusions are offered in Sec. VII. 



II. KINETIC MODELS OF LIQUID-CRYSTALLINE POLYMERS 



Consider a solution of liquid-crystalline polymers composed of rigid rodlike polymeric molecules. The kinetic models de- 
scribe the microstate of liquid-crystalline polymers by the one-particle distribution function ip(u; t), the probability density of 
finding a particle oriented along the unit vector u at time t. For simplicity, only spatially homogeneous systems are considered. 



A. Kinetic Equations of Polymer Dynamics 



Equations of motions for rigid rod models were proposed by Hess and Doi, flo, [/J |8p, and can be found in most of the textbooks 
on polymer kinetic theory (see e.g. [J3J [|, ||]). In the presence of a given homogeneous velocity gradient k, the time evolution of 
tp may be written as 

d t *jj = -i? u • [U X (K ■ UVO] + Ru ■ Dr4>Ru (l^rr) • (1) 

Here, R u = u x d/du is the rotational operator, d/du the gradient on the unit sphere, D r the rotational diffusivity and 6 /dip 
the Volterra functional derivative. The dimensionless free energy functional per molecule, F[ip] = Fo[ip] + F\ [tp], can be split 
into the entropy of molecular alignment, 

F [tp] = y"rfu^(u)ln^(u), (2) 
and the free energy contributions of interactions F\ . In the second virial approximation, F\ is given by 

FiM = £«/J(u,u')>>, (3) 

where v — n p d£ 2 denotes the reduced excluded volume of n p rods of length I and diameter d per unit volume. Here and below 
we use the following notations for averages: (/(u)) = Jdu/(u)?/'(u), ((/(u, u'))) = jjdudu'f(u, u')ip(u)ip(u'), where 
integration is performed over the three-dimensional unit sphere. 

The kinetic equation ([[]) is supplemented by the expression for the polymer contribution to the stress tensor [|[ ^o| , 

rP[^--3n p fc B T 1 ^[uxi? u ^-y]uy (4) 

In Eq. (Q), viscous contributions to t p are neglected since they are generally assumed to be negligible in liquid-crystalline 
polymers [||]. When Eq. (Q) is used in the Navier-Stokes equations, the latter, together with the kinetic equation ([l]), constitute 
a closed system describing the dynamics of the solution. 

The configuration dependent diffusion coefficient D r describes the hindrance of rotations due to neighboring rods. Following 
Doi and Edwards [ph, D r is approximated by 



l((Vl-(u-u') 2 



1 -2 



(5) 
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where D r , the rotational diffusion coefficient of a rod in an isotropic, semi-dilute solution of identical rods, is related to the 
rotational diffusion constant for a dilute solution, D r0 , by D T = cD r0 (n p £)~ 2 with an empirical coefficient c. It is generally 
believed that the self-consistent averaging approximation (Q) does not alter the dynamics ([[]) qualitatively [pT|. 

The nonlinearity of the kinetic equation ([l]) in tp reflects the mean-field character of the model that provides an effective 
one-body description of the many-rod system. The equilibrium distribution, ip CCi , corresponding to the stationary solution to the 
kinetic equation (jl|) in the absence of flow, n = 0, is given by tp cq oc exp [—[/]. The dimensionless self-consistent potential U 
is identified with 

U(u; i>) = = v J du'ftu, u')V-(u'), (6) 

where Eq. (^J) has been used. 

If only excluded-volume interactions are considered, the dimensionless second virial coefficient (3(u, u') for rigid rods was 
obtained by Onsager [55], 

/3(u, u') = |u x u'| = v/l-(u-u') 2 . (7) 

Eq. (0) has a simple geometric interpretation since |u x u'| is the area spanned by the vectors u and u'. The excluded-volume 
of two cylinders with length I, diameter d and orientations u and u' is therefore given by dl 2 1 u x u' | . Thus, F\ decreases as the 
polymers orient in the same direction and finally leads to the nematic phase when the effect of the excluded-volume becomes 
sufficiently strong. Note, that the second virial approximation becomes exact in the limit of high aspect ratio, Ijd — > oo. 

A variety of further expressions for F\ have been suggested in the literature [||. A particular simple, phenomenological 
expression was given by Maier and Saupe [ |23[ [J4j], 

/3 MS (u,u')=c -ci[(u.u') 2 -i], (8) 

which together with Eq. leads to 

F 1 MS ( r ^) = ^[c -c 1 , i^: aS 1 ], (9) 

where cq and c\ are parameters independent of tp. In Eq. (^), the important notion of the orientational order parameter a^T is 
introduced, = (uu — The isotropic state is characterized by = 0, while a^" ^ indicates (nematic) orientational 
ordering. 

The self-consistent equilibrium distribution corresponding to the potential (||) determines the isotropic nematic transition in 
equilibrium. The detailed form of the interaction potential can have significant effect on the behavior of the order parameter in 
the nematic phase. Specifically, the amount of order at the transition is known to be much smaller in the Maier-Saupe theory 
than in the Onsager model [pj[|. 

In the sequel, we consider generic interaction functionals of the form 

F 1 [^]=F 1 {M), (10) 



where M = {Mi, . . . , M mn } — { &2 , . . . , a.2 n } denotes the set of irreducible (anisotropic) moments of the distribution 
function up to order 2n, (a^ ) ai p 1 ...u n p Tl = (u ai u/s 1 . . . u an up n \, and B denotes the symmetric irreducible part of the 
tensor B [Eq]. In other words, the functional dependence of the interaction functional on the orientational distribution function 
ip comes only through the dependence of the moments of tp. For example, the tensor a^ contains five independent components. 
For convenience, these components are denoted by the scalars Mi, . . . , M5. The total number of components of irreducible 
moments up to order 2n is m n . In case n = 1 one has m„ = 5. In case n > 1, the components Mi for i > 5 contain the 
independent components of the tensors &2k , k = 2, . . . , n. 

The interaction functionals ( |l0| ) correspond to virial coefficients of the form /3(u, u') = /3((u • u') 2 — i) with polynomials 
/3(x) of degree n. The Maier-Saupe expression (^) is of the form ( |T(i| ) with n ~ 1. The Onsager potential (|7|), on the contrary, 
can only be cast into the form ( |lo| ) for n — > 00. In this case, systematic approximations of the form ( |Io| ) with finite n have been 



proposed in [EOJ. The lowest order approximation to the free energy functional (ph obtained in [ 20 1 for the Onsager potential is 



F 1 (- a 7) = -^^l-!- a 7: ai\ (11) 

Some evidence wasprovided in [ ^o| ] that the lowest order approximation represents a good approximation to F± in case of the 
Onsager potential (Q), at least on a representative subset of distribution functions. 
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Since the rotational diffusivity (^|) is related to the Onsager interaction functional F±, the same approximation can be employed 
to give [|§ 

D r = D[ n) (M). (12) 
Inserting the lowest order approximation for Fi one obtains [ p0[ | 

d[ 1] = (37r 2 /32)£) r [l - (3/2) a& : a?]" 1 . (13) 

The diffusion coefficient D, 1 ' is positive in the entire physically meaningful range of the order parameter . Expression ( |l3| ) 
should be compared with the Doi phenomenological result 

A-d = D t [l - (3/2) a^" : 'aT]- 2 . (14) 

For highly oriented nematics the expression ( pj[ ) is preferred compared to Doi's result [p7||. It should be stressed that our 
derivation of the diffusion coefficient does not need any further assumptions or adjustable parameters, while the derivation of 
Eq. (|l4|) requires the matching of D T , resp. F mt - MS in both, the isotropic and the fully ordered state [|| ^|. Since the rotational 
diffusivity D r is only determined within a phenomenological constant, the numerical factor 37r 2 /32 in Eq. ( fl3[ ) can safely be put 
equal to one. 



B. H-Theorem 



The time evolution of the free energy functional F = Fq + F\ is given by 



F =J d »W) dMu) - (15) 

Inserting the kinetic equation (jl|) into Eq. ([l5]) one obtains 

*_» ! ^_/ 4l "M. 1 U.S4y (16) 
\ Sip(u) Sip(u)/ 

In the absence of a velocity gradient, the dynamics ([I]) drives the system irreversibly to the unique equilibrium state ip oq with F 
being an ii-function of the relaxation dynamics. In the presence of a velocity gradient, Eq. (|6|) also describes the free energy 
exchange between the polymer and the solvent subsystem. The expression (Q) for the stress tensor can be rewritten in terms of 
partial rather then functional derivatives if the free energy F is of the form F = Fo + -Fi, with Fi given by Eq. ([l(]) depends on 
if> only through moments, 

t p = 3n p k B T , a^ + V] Au x R u (u ai Uj3 1 ■ ■ ■ u ak u,3 k ')}u\ - — . (17) 

k=1 X 1 0{ B. 2k ) ai l3 1 ---a k p k 

In the special case n = 1, Eq. ( |l7[ ) simplifies to 

f)F r)F 

r p = 3n p fc B T'a7 -2^^ -32 + 2^^: a 4 . (18) 
da 2 d& 2 

For the Maier-Saupe potential, Eq. ([)]), Eq. ( p^ ) coincides with the standard expression for the stress tensor given in [H]. 



III. CANONICAL DISTRIBUTION FUNCTIONS IN POLYMER DYNAMICS 



Within the kinetic models described in Sec. II, the dynamics and viscoelastic properties of liquid-crystalline polymers are 
given by Eqs. ([!]) and (Q). In viscoelastic flow calculations, the combined simulation of the Navier-Stokes equations coupled 
with the kinetic equation ([!]) is computationally very intensive. Therefore, there has been considerable interest in deriving 
sufficiently accurate closed form constitutive equations from kinetic theory. 
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A. The Closure Problem 

Isotropic, uniaxially and biaxially oriented phases of LCPs are conveniently described by the order parameter tensor . It 
is reasonable therefore to look for closed descriptions of the dynamics in terms of alone. The description of the viscoelastic 
properties of LCPs in terms of the stress tensor (|l7|), however, requires also the knowledge of moments up to order 2n + 2. The 
time evolution of moments a-y ' obtained from the kinetic equation ([[]) is of the form 

^i^7 =G<M(a£,... ,a^T) (19) 



with the functions 



G a?k~ajPj ( &2 ' ' ■ ' ' & 2i+2n ) = (^["«i H ft ' ' ' U aj U fij ] • (u X K • u) 

+ A (,l) (^['« Q1 «ft-s u ftl) 

+£)( ?l) ^ ( -Ru ['M ai ■ • ■ H aj M/3, ' ] • flu [u ai Ufa ■ ■ • Mfa ' ] ) T^^T ■ ( 20 ) 

Eqs. ( |l9| ) and (|o|) form a hierarchy of moment equations where moments of order 2j couple to moments of order (2j + 2n). 
Therefore, closed form constitutive equations in terms of low order moments cannot be derived exactly from the kinetic equation 
but require some approximations of the distribution function. There exist an enormous amount of closure approximations as 
well as many numerical tests of the approximations in various flow situations (see, e.g. || [30[ |T[ [32|]). 



B. Extremum Principle and Canonical Distribution Functions 

Following the approach described in [jl[], we here derive closed form constitutive equations for LCPs by applying the so-called 
quasi-equilibrium approximation. Within this approximation, canonical distribution functions are obtained from the extremum 
principle under constraints [Ol |l6|] 

F[i/>] -> min, 1= fduip(u), M k = J dum k (u)ip(u) (21) 

for k = 1, . . . , m n . Notice, that this principle does not coincide with the entropy maximum principle, because the free energy F, 
in general, is not proportional to the entropy. We consider in the sequel the special choice 'u ai u^ 1 ■ ■ ■ u ak up k for the functions 
mk(ii), such that M represents the set of the first n irreducible moments of the distribution function already introduced in 
Eq. (|To|). Inserting the functional form F = Fq + F\, where Fq is given by Eq. (||) and F\ by Eq. (|lo|), the solution to Eq. ( ^l| ) 
is given explicitly as 

m n 

r(u)=V eq (u)exp[^A fe m fc (u)+A ]. (22) 
fe=i 



The set of Lagrange multipliers A = {Ai, . . . , A m?i } = { A2 , . . . , A2„ } ensure the constraints in Eq. (|2l|), 

M k = / ciu m fc (u)f q (u) exp[^ A jJnj (u) + A ] , (23) 
J 3=1 

and Ao ensures the normalization of ip*. The interpretation of the CDF (^) as a quasi-equilibrium distributions has already 
been given in p3}|. 

The quasi-equilibrium free energy F* is defined as F* — F[ip*]. Inserting the form ([Z2|) one obtains 

F*(M) = J2^kM k + A . (24) 

k=i 



Note, that although the interaction part of the free energy is given by Eq. ( |10| ), the total quasi-equilibrium free energy F* cannot 
be expressed as a function of M_ explicitly. This difficulty results from the fact that the normalization Ao of the distribution ( p2| ) 
cannot be evaluated analytically in general. From Eq. (|24j), the Lagrange multipliers can be interpreted as the conjugate to the 
macroscopic variables, 

dF*(M) 

Afc = an/r — , k = 1, ...,m„. (25) 
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C. Macroscopic dynamics 

In the following, the canonical distribution functions ( p2| ) are employed in order to derive the macroscopic description of the 
polymer dynamics. We assume, that the CDF can be considered as representative states in the sense that a set of moments of 
the distribution function is approximated accuratelyby the corresponding moments evaluated with the CDF. A measure of the 
accuracy of this approximation is presented in Sec.jv] For improvements on this approximation see [ |l4] , [l7j|. 

The macroscopic time evolution is defined as the time evolution of the macroscopic variables evaluated with the CDF, M k = 
Jdu uik (u) dtip*. Inserting the kinetic equation ([j]) for dtip, one finds 



M fe = G fe (M) 



with 



(26) 



Gfe(M) = ( [R u m*(u)] • (u x k ■ u) 



where the symmetric, positive semi-definite matrix M is defined by 



M, 



A (,i) ([^um fe (u)].[i? u m,-(u)] 



(27) 



(28) 



Eqs. (P-6|), (|7|), ( P8| ) together with Eq. (22|) represent the closed set of macroscopic equations. The constitutive relation is 
obtained by evaluating the stress tensor m) with the CDF, t p '*(A1) = T p [ip*]. Here, it reads 



rP'*(M) - -3n p k B Tj2 ([u x R u m k (u)]u 



k=l 



dF*(M) 
a dM k 

The macroscopic free energy change, F*, is found from Eqs. ( |16[ ) and ( |26[ ) to be given by 

m n 
k,j=l 

and could have been obtained also by evaluating the free energy change ( |T6| ) on the CDF 
energy change (30) is negative semi-definite as is the underlying kinetic model. 



(29) 



(30) 



In the absence of flow, the free 



D. Alignment Tensor Models 

So far, the macroscopic description of LCP considered in the literature has almost exclusively considered the order parameter 
tensor as the only macroscopic variables. In our notation, this corresponds to the special choice of the order parameter tensor 
as the only macroscopic variable, n = 1 and = (uu). Due to symmetry, is the lowest non-trivial moment of ip. 

In this approximation, the dimensionless mean-field interaction potential derived from the free energy functional in case of 
the Onsager potential reads [Boh, 

m , . v 1 — |uu: a^ 

E/W(u, a 2 ) = -= 2 =, (31) 

v6*/l-|aj: a? 



which can be compared to the expression for the Maier-Saupe potential 

[/ms(u, a^ ) = c 2 - cii/uu: a^, (32) 

where c 2 is an arbitrary constant. 

In this approximation, the canonical distribution function (^) reduces to 

^a(u) =exp[-i//3'( r a7: + AT): uu + Aol/Z^ 1 ). (33) 

The Lagrange multipliers A 2 form a symmetric traceless matrix dual to a^ . The quasi-equilibrium distribution (B3) not 

. . * . . ss — 

only reduces to the equilibrium distribution for A 2 = but also gives the exact steady state solution A 2 = —j/2 for 
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homogeneous potential flows. Due to the occurrence of the second moment &2 > the quasi-equilibrium states ipX have to be 

determined self-consistently. Further manipulations are simplified by introducing new variables © = A2 — ^/3'( a^ : a^ ) 
and Go = Ao — In where © is a symmetric, traceless matrix. In terms of new variables, the quasi-equilibrium states ( f3~3| ) 
take the form 

ip* (u) = exp [u ■ ■ u + 9 ] , (34) 



where Go ensures the normalization of ip*. The distribution ([34|) is sometimes termed Bingham distribution. 

Accordingly, the diffusion coefficient D r is approximated by given in Eq. (|l3|). In the following, the macroscopic time 
evolution equations for a7 are obtained from Eqs. (26), (^7|), (||]) and read, 

a? =G* h (^) + G* d C^). (35) 
The presence of a flow field gives raise to the contribution G ; * in Eq. (|3~5|), 

G^Ca7) = n ■ + • k T + ^7 - 7: a* 4 , (36) 

where 7 = k + k t denotes the rate-of-strain tensor and a^ = (uuuu) A is evaluated with the Bingham distribution (j34}). The 
contribution G*j of the Brownian motion and the interaction potential is given by 

G^) = M*:A7, (37) 

where A2 = dF* (a^ ) /d'a^ and the matrix M* is given by 



In Eq. (38), use has been made of the fact that the approximated diffusion coefficient D r is independent of the orientation u. 
Inserting the expression for the free energy one arrives at 

G$(ar) =-6Df } a7 + 6$ 1) v0'(&: a?) (a? • a? + ~"a7 - a? : aj), (38) 

The Bingham distribution ([34|) has already been employed in the literature [ |30[ f34| ] to derive closed form constitutive equations 
from the kinetic model (jl|). However, the expression for G^ differs from the results given in [[30, 34 1 not only in the diffusion 

coefficient D^ 1 . It also generalizes earlier result to general interaction potentials of the form (|To|). For the approximation ( f3l| ) to 
the Onsager potential, for example, a non-polynomial dependence on the order parameter a^ occurs, which becomes important 
in the nematic state |20|. Comparative studies of various closure approximations to the kinetic model for the Maier-Saupe 
potential have clearly demonstrated that closures based on the Bingham distribution are preferable compared to more ad hoc 
closure approximations [^8j] . 

In the present case, the macroscopic free energy ( p4| ) reduces to 

F*(a£) = A7: a7 + A . (39) 
The macroscopic polymer contribution to the stress tensor ( p^ ) simplifies to 

tP<* = 2k B T(^ • AT -a*: A^) . (40) 

Expression fiQ) is also obtained when the macroscopic entropy (p9[) is inserted into Eq. (Eh. 



IV. NUMERICAL INTEGRATION SCHEME 

Eqs. d26| ) are the closed form macroscopic equations. Together with the expression (^) for the stress tensor they provide 
the macroscopic constitutive equation of the system. However, the time evolution equations ( |26| ) still contain the Lagrange 
multipliers, that have to be determined such that the constraints in Eq. ( j2l| ) are satisfied. 

To deal with this situation, a numerical integration scheme for the macroscopic equations in the context of dilute solutions of 
flexible polymers has been proposed in Ref. [El], This scheme is general enough to be applied also to the present situation. We 
present here this integration scheme only briefly, referring the reader to ^ for further details. 
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A. Dynamics of Dual Variables 



In the extremum principle (|2JJ), the Lagrange multipliers occur naturally to satisfy the constraints. Due to Eq. (|25|), the 
Lagrange multipliers are interpreted as conjugate variables. Instead of eliminating the Lagrange multipliers A in favor of the 
macroscopic variables A = A(M), the Lagrange multipliers play the role as dynamical variables that determine the values of 



the macroscopic variables, M = A/ (A). Note, that the functions M(A) are given by Eq. (23). With the help of the Legendre 



transform of the macroscopic free energy, G(A) = F* (M_) + J2k=i Afc-M/., the time evolution of A is given by 

Afe = E( C )*7 A ^A), (41) 
where C _1 is the inverse of the matrix C k j = <9 2 G(A) / dA k dAj and Mj (A) is given by the right hand side of Eq. (|26|). 



B. Numerical Integration Scheme 



The reformulation of the macroscopic dynamics ( p6[ ) in terms of the dual variables A described in Sec. IV A is suitable for 
numerical implementation. The Lagrange multipliers A now play the role of independent dynamic variables, instead of M. In 
order to advance given values A(t) at time t to their values A(t + r), with small time step r, the following first order integration 
scheme is proposed in Ref . [[[J : 

1 . The new values of the macroscopic variables M_(t + t) are found from the values A(t) by 

M fc (t + T)-Af fc (t) (42) 
r 

where G k ( M ) denotes the right hand side of Eq. ( p6[ ) with A(t) the actual values of the Lagrange multipliers. 

2. The matrix Ckj is evaluated from 

C kj (t) = {m k m 3 ) K (t) - M k {t)M.j(t), (43) 
where (m k mi) A is calculated with the distribution function ( |22| ) where A = A{t). 

3. The n x n matrix C(<) is inverted numerically to give C _1 (i). 

4. The values of the Lagrange multipliers, A(t + r), are given by 

A k (t + r)-A k (t) = g c _ x (t) Mjjt + ^-Mjjt) (44) 



This concludes one time step of integration. The integration scheme has to be supplemented by initial conditions A(0). The 
special case of equilibrium initial conditions corresponds to A(0) = 0. The appearance of the correlation matrix C as a tool 
to recompute moments into Lagrange multipliers in this scheme is intimately related to the structure of the quasi-equilibrium 
approximation. 

Note, that evaluating the matrix C requires moments of ^ which are of higher order than the macroscopic variables them- 
selves. If n denotes the number of macroscopic variables, the numerical integration scheme requires in total n(n + 3)/2 integrals 
per time step to evaluate M and C. Due to the symmetry of C, this number is of order n 2 /2 rather than n 2 . The evaluation of 
all these integrals with standard numerical methods might be time-consuming, especially for high-dimensional integrals. It is 
demonstrated in [[35]] that under certain circumstances these integrals can be evaluated efficiently by adapting methods of fast 
Fourier transformations. This computational issue is out of scope of this paper. 



C. Numerical Implementation of Alignment Tensor Models 



For the special case of alignment tensor models described in Sec. [ED, the general integration scheme outlined in Sec. IV B 



is applied. The integrals over the unit sphere associated with the evaluation of moments of the canonical distribution function 
are performed numerically. Different resolutions of the numerical integration are used to obtain more accurate results by ex- 
trapolating to infinite resolution. The normalization of u, u 2 = 1, and the consistency relations (uuuu) A : 1 = (uu) A are 
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used as additional checks of the results of the numerical integration. The inversion of the matrix C is done by using the explicit 
expression for the inverse of a symmetric 6x6 matrix as given by the symbolic computation software Mathematica™. The 
Lagrange multipliers are advanced by Eq. (|4|). The resulting algorithm has been implemented and tested under various circum- 
stances. This numerical implementation differs from the one proposed in [p0||. In the latter reference, approximate expressions 
for the eigenvalues of &% as a function of the eigenvalues of are proposed, whereas in the present approach the necessity of 
expressing the Lagrange multipliers A in terms of the moments is circumvented by the scheme (|4^)-(p4^. Together with suitable 
rotations to diagonalizing coordinate systems, these expressions allow direct numerical integration of Eq. (B5b. 



V. ACCURACY OF THE APPROXIMATION BY CDF 



The use of the CDFs for the macroscopic description described in Sec. Ill C imposes approximations to the underlying dy- 
namics. Therefore, the accuracy of this approximation needs to be discussed. 

The accuracy of the approximation by CDFs can be studied in two ways: First, the result of integrating the macroscopic 
dynamics can be compared to direct numerical integration (e.g. by Brownian dynamics) of the underlying kinetic model. Alter- 
natively, the dynamic variance A that plays a central role in the method of invariant manifold [|l4| may be used in order to access 
the accuracy of the approximation. The advantage of using A to estimate the accuracy of the approximation is that it can be 
evaluated while integrating the macroscopic dynamics, without performing numerics for the underlying microscopic dynamics. 
For the use of the dynamic variance in a mixed micro-macro computation see []l7p. 

The dynamic variance A is defined as 

A(u;M) = Jf-^M;, (45) 

3=1 3 

where Jip* denotes the right hand side of Eq. ([!]) evaluated with the CDF ip*, 

Jr = -R u ■ [u x (#6- w/>*)] + D r Y,Ru ■ ip*R u m k (u)A k . (46) 

k 

The dynamic variance A can be interpreted as the difference of the microscopic time evolution evaluated with ip*, Eq. (^6|), and 
the time evolution of the CDFs due to the macroscopic dynamics. If A = 0, the CDFs are said to form an invariant manifold 
of the dynamics ([!]). In case A / 0, the method of invariant manifold uses A to obtain improved manifolds of distribution 
functions. Here and in a companion paper [JTf], | A| is used as a measure of the accuracy of the approximation by CDFs, where 
the norm | • | is specified later. 

Consider functions Yk(M) defined by Yk(M) = Jduyk(u)A(u;M_). Some choices of yk(u) are discussed below. The 
functions Y k can be interpreted as the difference of the time evolution of the macroscopic moments if evaluated from the 
micr oscop ic dynamics and from the macroscopic dynamics. For the special case of alignment tensor models, introduced in 
Sec. Ill D| , functions Y k found from Eq. ( p3| ) read 



Y k = (y fe (u)uu: [3k - 2© ■ P • k + 4D r • P • A 2 - 6D r A 2 

\ '/A 

- fe(u)(uu)) A -{ft(u)) A (uu} A ]: -— : a 2 , (47) 

a a 2 

where P = 1 — uu is the projector perpendicular to u, is given by Eq. ([?5|) and 

= CT 1 - ufi'll - 2v(3"^ a? . (48) 

<9a 2 

We have Yo = for yo = 1, since the first and second term in Eq. (^) conserve the normalization of the distribution function 
independently. Further, Y k vanishes by construction if yk corresponds to components of the tensor uu. For y k corresponding 
to uuuu or higher order tensors, the quantities Y k contain valuable information about the variance A and may be used as a 
measure for I A I . 



VI. SOME NUMERICAL RESULTS 



In the present section, the macroscopic equations ( |37[ ) with ( [34] ) are integrated numerically using the dual integrator presented 
in Sec. tVBl. 
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A. Some Numerical Results for Steady Shear Flow 

The integration scheme described above was used to obtain the time evolution of the macroscopic variables from Eq. (B7| 



in steady shear flow. The functions G and F are given in Eqs. (36) and (f38|), respectively with n a p = A /5 a i5p2- We employ 



the approximation (31 ) to the Onsager excluded volume potential and expression ( |13| ) for the diffusion coefficient. The minimal 
value V2 of the nematic potential that allows only a stable nematic phase was estimated from steady state values of the numerical 
solution of Eq. ([37) in the absence of flow, G = 0. For almost isotropic initial conditions a value v 2 ~ 6.5 was obtained. 



Following Larson [36|, we set v 2 equal to this minimal value. Fig. [j] shows the dimensionless shear viscosity as a function of 
shear strain for two low values of the dimensionless shear rate Pe as obtained form the numerical integration for equilibrium 
initial condition. The Pec let number Pe is defined as the dimensionless shear rate Pe = j/6D*, where D* is the diffusion 
coefficient D T defined in Eq. (fy and below at the concentration v%. It is obvious from Fig. [j] that the present macroscopic 
description shows a so-called 'tumbling' or 'wagging' regime [[36]]. 

Rheological and rheo-optical experiments on lyotropic solutions of rodlike polymers have clearly demonstrated the existence 
of a 'tumbling' or 'wagging' regime for low shear rates. In Fig. pi th e experimental result for the shear viscosity of the lyotropic 



nematic polymer n-octylcyanobiphenyl is given as reported in [p7|]. Note, that the experimental results of Refs. [|37| |3gj are 
obtained from birefringence data of a monodomain liquid-crystalline polymer. Therefore, neglecting spatial inhomogeneities in 
the kinetic model is justified. The experimental data were obtained in a parallel-plate flow from birefringence measurements 
for a shear rate of 7 = 16 s -1 . Qualitatively, the simulation results show similar behavior than the experimental findings. 
However, the damping of the oscillations of the shear viscosity is not seen in Fig. [j]. The authors of Ref. J3~^ ] speculate that 
this behavior might be due to out-of-plane director orientations. Another dissipative mechanism not included in the present 
kinetic model is translational diffusion which could also explain the damping of the oscillations in Fig. [I]. We mention, that 
the authors of Ref. [|8j] obtained a fit to experimental data on a similar liquid-crystalline polymer which is comparable to our 
simulation results. Their procedure requires fitting of three parameters of a generalized version of Ericksen's theory of nematic 
liquid crystals | 



VII. CONCLUSION 



A systematic approach to the derivation and implementation of constitutive equations proposed in |Jl|] is applied in the present 
work to the dynamics of liquid-crystalline polymers. Employing the quasi-equilibrium approximation, a set of canonical distri- 
bution functions is obtained which is further used to derive constitutive equations. The numerical implementation and a measure 
of the accuracy of the approximation is discussed. The present approach is illustrated for the kinetic model of liquid-crystalline 
polymers with the Onsager excluded volume potential in steady shear and steady elongational flow. 
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FIG. 1: Time evolution of the dimensionless shear viscosity as a function of shear strain. The nematic potential was chosen to be vi — 6.5. 
Solid and dashed line correspond to the dimensionless shear rate Pe = 5 and Pe = 10, respectively. 
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FIG. 2: Time evolution of the shear viscosity of n-octylcyanobiphenyl as a function of shear strain for shear rate j — 16 s 1 . Data are taken 
from[0. 
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